Fast followers?/Speed contagion: assessing the impact of Montreal F1 Grand Prix on high-speed ticketing rates (2000-2022)

Step 3. Check analytical strategies (part 2)

Author

Andrés González Santa Cruz

Published

July 7, 2025

Code
# remove objects and memory
rm(list=ls());gc()
          used (Mb) gc trigger (Mb) max used (Mb)
Ncells  835948 44.7    1627021 86.9  1178085   63
Vcells 1686107 12.9    8388608 64.0  3406660   26
Code
#remove images
while(!dev.cur())dev.off()
cat("\014")
Code
load(paste0(getwd(),"/_data/step2.RData"))

Load libraries and data

Code
#borrar caché
#system("fc-cache -f -v")

#check R version
if(Sys.info()["sysname"]=="Windows"){
if (getRversion() != "4.4.1") { stop("Requiere versión de R 4.4.1. Actual: ", getRversion()) }
}
if(Sys.info()["sysname"]=="Linux"){
if (getRversion() != "4.4.1") { stop("Requiere versión de R 4.4.1. Actual: ", getRversion()) }
}
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:

# install.packages(c("dplyr", #for data
#                    "tidyr", #for data
#                    "lubridate",  #for dates
#                    "openxlsx", #for excel files
#                    "rio", #for importing and exporting data
#                    "purrr", #for iterating in databases
#                    "devtools", #for external packages
#                    "DiagrammeR", #para visualizar DAG
#                    "dagitty",
#                    "ggdag",
#                    "ggplot2" #for graphics
#                    "kableExtra", #pretty tables
#                    "quarto", #for documents
#                    "geosphere" #for coordinates and classifying
#                    "geepack", #for regression
#                    "glmmTMB", #For GLMMs
#                    "DHARMa", #For residual diagnostics
#                    "car", #For hypothesis testing
#                    "brms", #Bayesian model
#                    "bayesplot",
#                    "loo",
#                    "Synth", #for synthetic control method
#                    "weathercan",#for weather data
#                    "sandwich", #cluster robust intervals
#                    "emmeans", #for predictions
#                    "gnm", #Conditional Poisson models
#                    "splines", #nonlinearity 
#                    "geeM", #negative binomial and more flexible GEE models
#                    "PanelMatch", #Matching technique with panel data
#                    "scpi" #control sintético
#                    "nixtlar", #for time series analysis and prediction
#                    "CausalImpact" #for time series causal impact
#                    "forecast" #for time series analysis prediction and decomposition
#                    ))
library(dplyr); library(lubridate); library(tidyverse); library(openxlsx); library(rio); library(purrr); library(dagitty); library(ggdag); library(kableExtra); library(geosphere); library(geepack); library(lme4); library(glmmTMB); library(DHARMa); library(car); library(brms); library(bayesplot); library(loo); library(Synth);  library(weathercan); library(sandwich); library(emmeans); library(gnm); library(splines); library(geeM); library(plm); library(PanelMatch); library(scpi); library(fect); library(nixtlar); library(CausalImpact); library(forecast)

Adjuntando el paquete: 'dplyr'
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

Adjuntando el paquete: 'lubridate'
The following objects are masked from 'package:base':

    date, intersect, setdiff, union
Warning: package 'readr' was built under R version 4.4.3
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ forcats 1.0.0     ✔ stringr 1.5.1
✔ ggplot2 3.5.2     ✔ tibble  3.2.1
✔ purrr   1.0.4     ✔ tidyr   1.3.1
✔ readr   2.1.5     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Adjuntando el paquete: 'dagitty'


The following object is masked from 'package:rio':

    convert



Adjuntando el paquete: 'ggdag'


The following object is masked from 'package:stats':

    filter



Adjuntando el paquete: 'kableExtra'


The following object is masked from 'package:dplyr':

    group_rows


Cargando paquete requerido: Matrix


Adjuntando el paquete: 'Matrix'


The following objects are masked from 'package:tidyr':

    expand, pack, unpack



Adjuntando el paquete: 'lme4'


The following object is masked from 'package:rio':

    factorize


This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')

Cargando paquete requerido: carData


Adjuntando el paquete: 'car'


The following object is masked from 'package:purrr':

    some


The following object is masked from 'package:dplyr':

    recode


Cargando paquete requerido: Rcpp

Loading 'brms' package (version 2.22.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').


Adjuntando el paquete: 'brms'


The following object is masked from 'package:glmmTMB':

    lognormal


The following object is masked from 'package:lme4':

    ngrps


The following object is masked from 'package:stats':

    ar


This is bayesplot version 1.12.0

- Online documentation and vignettes at mc-stan.org/bayesplot

- bayesplot theme set to bayesplot::theme_default()

   * Does _not_ affect other ggplot2 plots

   * See ?bayesplot_theme_set for details on theme setting


Adjuntando el paquete: 'bayesplot'


The following object is masked from 'package:brms':

    rhat


This is loo version 2.8.0

- Online documentation and vignettes at mc-stan.org/loo

- As of v2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session. 

- Windows 10 users: loo may be very slow if 'mc.cores' is set in your .Rprofile file (see https://github.com/stan-dev/loo/issues/94).

##
## Synth Package: Implements Synthetic Control Methods.


## See https://web.stanford.edu/~jhain/synthpage.html for additional information.
Warning: package 'weathercan' was built under R version 4.4.3
As of v0.7.2, the `normals` column in `stations()` reflects whether or not there
are *any* normals available (not just the most recent).
Welcome to emmeans.
Caution: You lose important information if you filter this package's results.
See '? untidy'
Registered S3 method overwritten by 'lfe':
  method    from 
  nobs.felm broom

Adjuntando el paquete: 'plm'

The following objects are masked from 'package:dplyr':

    between, lag, lead


Adjuntando el paquete: 'PanelMatch'

The following object is masked from 'package:tidyr':

    extract

The following object is masked from 'package:stats':

    weights

Registered S3 method overwritten by 'GGally':
  method from   
  +.gg   ggplot2
Cargando paquete requerido: bsts
Cargando paquete requerido: BoomSpikeSlab
Cargando paquete requerido: Boom

Adjuntando el paquete: 'Boom'

The following objects are masked from 'package:brms':

    ddirichlet, rdirichlet

The following object is masked from 'package:stats':

    rWishart


Adjuntando el paquete: 'BoomSpikeSlab'

The following object is masked from 'package:stats':

    knots

Cargando paquete requerido: zoo

Adjuntando el paquete: 'zoo'

The following objects are masked from 'package:base':

    as.Date, as.Date.numeric

Cargando paquete requerido: xts

######################### Warning from 'xts' package ##########################
#                                                                             #
# The dplyr lag() function breaks how base R's lag() function is supposed to  #
# work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
# source() into this session won't work correctly.                            #
#                                                                             #
# Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
# conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
# dplyr from breaking base R's lag() function.                                #
#                                                                             #
# Code in packages is not affected. It's protected by R's namespace mechanism #
# Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
#                                                                             #
###############################################################################

Adjuntando el paquete: 'xts'

The following objects are masked from 'package:dplyr':

    first, last


Adjuntando el paquete: 'bsts'

The following object is masked from 'package:BoomSpikeSlab':

    SuggestBurn

Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

Adjuntando el paquete: 'forecast'

The following object is masked from 'package:brms':

    ma
Code
#special repository indicated or the package
if(!require(weathercan)){
   install.packages("weathercan", 
                  repos = c("https://ropensci.r-universe.dev", "https://cloud.r-project.org")); library(weathercan)
  }

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

if(!require(bpmn)){devtools::install_github("bergant/bpmn")}
Cargando paquete requerido: bpmn
Code
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
sum_dates <- function(x){
 
  cbind.data.frame(
    min= as.Date(min(unclass(as.Date(x)), na.rm=T), origin = "1970-01-01"),
    p001= as.Date(quantile(unclass(as.Date(x)), .001, na.rm=T), origin = "1970-01-01"),
    p005= as.Date(quantile(unclass(as.Date(x)), .005, na.rm=T), origin = "1970-01-01"),
    p025= as.Date(quantile(unclass(as.Date(x)), .025, na.rm=T), origin = "1970-01-01"),
    p25= as.Date(quantile(unclass(as.Date(x)), .25, na.rm=T), origin = "1970-01-01"),
    p50= as.Date(quantile(unclass(as.Date(x)), .5, na.rm=T), origin = "1970-01-01"),
    p75= as.Date(quantile(unclass(as.Date(x)), .75, na.rm=T), origin = "1970-01-01"),
    p975= as.Date(quantile(unclass(as.Date(x)), .975, na.rm=T), origin = "1970-01-01"),
    p995= as.Date(quantile(unclass(as.Date(x)), .995, na.rm=T), origin = "1970-01-01"),
    p999= as.Date(quantile(unclass(as.Date(x)), .999, na.rm=T), origin = "1970-01-01"),
    max= as.Date(max(unclass(as.Date(x)), na.rm=T), origin = "1970-01-01")
  )
}
smd_bin <- function(x,y){
  z <- x*(1-x)
  t <- y*(1-y)
  k <- sum(z,t)
  l <- k/2
  
  return((x-y)/sqrt(l))
  
}

theme_custom_sjplot2 <- function(base_size = 12, base_family = "") {
  theme_minimal(base_size = base_size, base_family = base_family) +
    theme(
      # Text elements
      text = element_text(size = base_size, family = base_family),
      plot.title = element_text(face = "bold", hjust = 0.5, size = base_size * 1.2),
      plot.subtitle = element_text(hjust = 0.5, margin = margin(b = 10)),
      axis.title = element_text(size = base_size, face = "bold"),
      axis.text = element_text(size = base_size * 0.8),
      axis.text.x = element_text(angle = 0, hjust = 0.5, vjust = 0.5),
      axis.text.y = element_text(angle = 0, hjust = 1, vjust = 0.5),
      axis.title.x = element_text(margin = margin(t = 10)),
      axis.title.y = element_text(margin = margin(r = 10)),
      
      # Plot layout
      plot.margin = margin(t = 20, r = 20, b = 20, l = 20),
      panel.grid.major = element_line(color = "grey80"),
      panel.grid.minor = element_blank(),
      legend.position = "right",
      legend.text = element_text(size = base_size * 0.8),
      legend.title = element_text(size = base_size, face = "bold"),
      legend.background = element_rect(fill = "white", colour = NA),
      legend.box.background = element_rect(colour = "grey80", linetype = "solid"),
      legend.key = element_rect(fill = "white", colour = "white")
    )
}

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
num_cores <- parallel::detectCores() -1
data.table::setDTthreads(threads = num_cores)#restore_after_fork = NULL, throttle = NULL)

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#CONFIG #######################################################################
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

options(scipen=2) #display numbers rather scientific number


nixtlar::nixtla_client_setup(api_key = readLines(paste0(gsub("f1/","", getwd() ),"/key.txt"))[[2]])
Warning in readLines(paste0(gsub("f1/", "", getwd()), "/key.txt")): incomplete
final line found on 'H:/Mi unidad/PERSONAL
ANDRES/UCH_salud_publica/pasantia/f1/key.txt'
API key has been set for the current session.

Forecast with Time GPT

Code
#https://cran.r-project.org/web/packages/nixtlar/vignettes/prediction-intervals.html
#https://cran.r-project.org/web/packages/nixtlar/vignettes/exogenous-variables.html
#https://cran.r-project.org/web/packages/nixtlar/vignettes/cross-validation.html
#https://cran.r-project.org/web/packages/nixtlar/vignettes/anomaly-detection.html

collisions_weather_corr_rect_synth_tr$Y_lag10 <- 
collisions_weather_corr_rect_synth_tr|> 
    mutate(Y_lag10 = dplyr::lag(rate_veh, 10, order_by = yday_corr))|>
    mutate(Y_lag10 = dplyr::lag(rate_veh, 10, order_by = yday_corr))|> 
    ungroup()|>
    pull(Y_lag10) 

collisions_weather_corr_rect_synth_tr$Y_lag10 <- 
collisions_weather_corr_rect_synth_tr|>
    mutate(
        third_value = nth(na.omit(Y_lag10), 1),
        Y_lag10 = ifelse(is.na(Y_lag10), third_value, Y_lag10)
    )|>
    dplyr::select(Y_lag10)|>
    ungroup()|>
    pull(Y_lag10)


df35 <- collisions_weather_corr_rect_synth_tr %>% 
  arrange(tr_contr_sens, year.x, yday_corr) %>%                # 1,2,3… days relative to race
  filter(yday_corr<36) |> 
  transmute(
    unique_id = "treated",             # exactly this name
    ds        = as.POSIXct(as.Date("2019-05-19") + (yday_corr - 1)),
    y         = rate_veh,
    mean_median_lag_2_prec_median_imp         = mean_median_lag_2_prec_median_imp,
    mean_median_total_precip_median_lin= mean_median_total_precip_median_lin,
    mean_min_temp_mean_lin=  mean_min_temp_mean_lin,
    mean_max_temp_mean_lin= mean_max_temp_mean_lin,
    Y_lag10= Y_lag10
  )
df35_post <- collisions_weather_corr_rect_synth_tr %>% 
  arrange(tr_contr_sens, year.x, yday_corr) %>%                # 1,2,3… days relative to race
  filter(yday_corr>=36, yday_corr<=max_time) |>                     # take the first 30
  transmute(
    unique_id = "treated",             # exactly this name
    ds        = as.POSIXct(as.Date("2019-05-19") + (yday_corr - 1)),
    y         = rate_veh,
    mean_median_lag_2_prec_median_imp         = mean_median_lag_2_prec_median_imp,
    mean_median_total_precip_median_lin= mean_median_total_precip_median_lin,
    mean_min_temp_mean_lin=  mean_min_temp_mean_lin,
    mean_max_temp_mean_lin= mean_max_temp_mean_lin,
    Y_lag10= Y_lag10
  )


df33 <- collisions_weather_corr_rect_synth_tr %>% 
  arrange(tr_contr_sens, year.x, yday_corr) %>%                # 1,2,3… days relative to race
  filter(yday_corr<34) |> 
  transmute(
    unique_id = "treated",             # exactly this name
    ds        = as.POSIXct(as.Date("2019-05-19") + (yday_corr - 1)),
    y         = rate_veh,
    mean_median_lag_2_prec_median_imp         = mean_median_lag_2_prec_median_imp,
    mean_median_total_precip_median_lin= mean_median_total_precip_median_lin,
    mean_min_temp_mean_lin=  mean_min_temp_mean_lin,
    mean_max_temp_mean_lin= mean_max_temp_mean_lin,
    Y_lag10= Y_lag10
  )
df33_post <- collisions_weather_corr_rect_synth_tr %>% 
  arrange(tr_contr_sens, year.x, yday_corr) %>%                # 1,2,3… days relative to race
  filter(yday_corr>=34, yday_corr<=max_time) |>                     # take the first 30
  transmute(
    unique_id = "treated",             # exactly this name
    ds        = as.POSIXct(as.Date("2019-05-19") + (yday_corr - 1)),
    y         = rate_veh,
    mean_median_lag_2_prec_median_imp         = mean_median_lag_2_prec_median_imp,
    mean_median_total_precip_median_lin= mean_median_total_precip_median_lin,
    mean_min_temp_mean_lin=  mean_min_temp_mean_lin,
    mean_max_temp_mean_lin= mean_max_temp_mean_lin,
    Y_lag10= Y_lag10
  )



fc <- nixtla_client_forecast(
  df35,
  h     = 7,
  level = c(95),
  freq  = "D",          # <— this fixes the infer_frequency error
  finetune_steps = 10, #Number of steps used to finetune 'TimeGPT' in the new data.
  finetune_loss ="rmse"
)
Using historical exogenous features: mean_median_lag_2_prec_median_imp, mean_median_total_precip_median_lin, mean_min_temp_mean_lin, mean_max_temp_mean_lin, Y_lag10
Code
fc33 <- nixtla_client_forecast(
  rbind(mutate(df33[1:2,], ds=ds-3), df33),
  h     = 7,
  level = c(95),
  freq  = "D",          # <— this fixes the infer_frequency error
  finetune_steps = 10, #Number of steps used to finetune 'TimeGPT' in the new data.
  finetune_loss ="rmse"
)
Using historical exogenous features: mean_median_lag_2_prec_median_imp, mean_median_total_precip_median_lin, mean_min_temp_mean_lin, mean_max_temp_mean_lin, Y_lag10
Code
p <- nixtla_client_plot(rbind(df35, df35_post), fc,max_insample_length = 100)
p33 <- nixtla_client_plot(rbind(mutate(df33[1:2,], ds=ds-3), df33, df33_post), fc33,
                        max_insample_length = 100)


# 1) Banda de CI (layer 1)
p$layers[[1]]$geom_params$fill  <- "#FFC10740"  # amarillo semitransparente
p$layers[[1]]$geom_params$alpha <- 0.4

# 2) Línea histórica y pronóstico (layer 2)
# p$layers[[2]]$aes_params$colour  <- "steelblue"
p$layers[[2]]$aes_params$linewidth <- 1.0


p33$layers[[1]]$aes_params$fill  <- "#eb8f8f"  # verde semitransparente
p33$layers[[1]]$aes_params$alpha <- 0.5


# 1) Banda de CI (layer 1)
p33$layers[[1]]$geom_params$fill  <- "#FFC10740"  # amarillo semitransparente
p33$layers[[1]]$geom_params$alpha <- 0.4

# 2) Línea histórica y pronóstico (layer 2)
# p$layers[[2]]$aes_params$colour  <- "steelblue"
p33$layers[[2]]$aes_params$linewidth <- 1.0


p33$layers[[1]]$aes_params$fill  <- "#eb8f8f"  # verde semitransparente
p33$layers[[1]]$aes_params$alpha <- 0.5

t0 <- min(p$data$ds)
t0_33 <- min(p33$data$ds)
t0_posix <- min(p$data$ds)
t0_33_posix <- min(p33$data$ds)

xint35   <- t0_posix + days(34)
t0 <- as.Date(min(p$data$ds))
t0_33 <- as.Date(min(p33$data$ds))

# 3) Aplica scale_x_datetime() para transformar POSIXct → numérico
p +
  scale_color_manual(
    values = c(
     # y        = "#FFD600",   # amarillo intenso
      TimeGPT  = "red"    # violeta (DarkOrchid / BlueViolet)
    ), 
    labels= c( TimeGPT = "Prediction (95% IC in gray area)" )
  ) +
    scale_x_datetime(
        name   = "Days of follow-up (35 days before race,\n7 days after the race)",
        limits = c(min(p$data$ds), max(p$data$ds)), 
        breaks = seq(min(p$data$ds), max(p$data$ds), by = "7 days"),
        labels = function(x) {
            # convierto POSIXct a Date y resto t0
            as.integer(as.Date(x) - t0 + 1)
        }
    ) +
    theme_classic(base_size = 13) +
    theme(
        strip.text   = element_blank(),   # quita el label de faceta
        strip.background = element_blank(),
        legend.position= "bottom"
    )+
      geom_vline(
    xintercept = xint35,
    linetype   = "dashed",
    color      = "black",
    linewidth       = 1
  ) +
  labs(y="Highspeed collisions per 1M vehicles")
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Prediciton on the treated (adjusting on min max temperature and lagged precipitiation (2 days) and lagged rates

Prediciton on the treated (adjusting on min max temperature and lagged precipitiation (2 days) and lagged rates

Code
p33+
    scale_color_manual(
        values = c(
            # y        = "#FFD600",   # amarillo intenso
            TimeGPT  = "red"    # violeta (DarkOrchid / BlueViolet)
        ), 
        labels= c( TimeGPT = "Prediction (95% IC in gray area)" )
    ) +
    scale_x_datetime(
        name   = "Days of follow-up (33 days before race,\n7 days after the race)",
        limits = c(min(p33$data$ds), max(p33$data$ds)), 
        breaks = seq(min(p33$data$ds), max(p33$data$ds), by = "7 days"),
        labels = function(x) {
            # convierto POSIXct a Date y resto t0
            as.integer(as.Date(x) - t0_33 + 1)
        }
    ) +
    theme_classic(base_size = 13) +
    theme(
        strip.text   = element_blank(),   # quita el label de faceta
        strip.background = element_blank(),
        legend.position= "bottom"
    )+
    geom_vline(
        xintercept = xint35-days(2),
        linetype   = "dashed",
        color      = "black",
        linewidth       = 1
    ) +
    labs(y="Highspeed collisions per 1M vehicles")
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Prediciton on the treated (adjusting on min max temperature and lagged precipitiation (2 days) and lagged rates

Prediciton on the treated (adjusting on min max temperature and lagged precipitiation (2 days) and lagged rates

Only Quebec as controls

Code
## ── Lags ──
collisions_weather_corr_rect_quebec_synth_tr$Y_lag10 <-
  collisions_weather_corr_rect_quebec_synth_tr |>
    arrange(tr_contr, year.x, yday_corr) |>
    mutate(
      Y_lag10 = dplyr::lag(rate_veh, 10, order_by = yday_corr),
      Y_lag10 = ifelse(is.na(Y_lag10),
                       first(na.omit(Y_lag10)),
                       Y_lag10)
    ) |>
    ungroup() |>
    pull(Y_lag10)

## ── Data frames (35 y 35 + 7) ──
df35_qc <- collisions_weather_corr_rect_quebec_synth_tr %>% 
  arrange(tr_contr, year.x, yday_corr) %>%
  filter(yday_corr < 36) %>%
  transmute(
    unique_id = "treated",
    ds        = as.POSIXct(as.Date("2019-05-19") + (yday_corr - 1)),
    y         = rate_veh,
    mean_median_lag_2_prec_median_imp,
    mean_median_total_precip_median_lin,
    mean_min_temp_mean_lin,
    mean_max_temp_mean_lin,
    Y_lag10
  )

df35_post_qc <- collisions_weather_corr_rect_quebec_synth_tr %>% 
  arrange(tr_contr, year.x, yday_corr) %>%
  filter(yday_corr >= 36, yday_corr <= max_time) %>%
  transmute(
    unique_id = "treated",
    ds        = as.POSIXct(as.Date("2019-05-19") + (yday_corr - 1)),
    y         = rate_veh,
    mean_median_lag_2_prec_median_imp,
    mean_median_total_precip_median_lin,
    mean_min_temp_mean_lin,
    mean_max_temp_mean_lin,
    Y_lag10
  )

## ── TimeGPT ──
fc_qc <- nixtla_client_forecast(
  df35_qc,
  h     = 7,
  level = c(95),
  freq  = "D",
  finetune_steps = 10,
  finetune_loss  = "rmse"
)
Using historical exogenous features: mean_median_lag_2_prec_median_imp, mean_median_total_precip_median_lin, mean_min_temp_mean_lin, mean_max_temp_mean_lin, Y_lag10
Code
p_qc <- nixtla_client_plot(
  rbind(df35_qc, df35_post_qc),
  fc_qc,
  max_insample_length = 100
)

## Personalización idéntica a tu original
p_qc$layers[[1]]$geom_params$fill   <- "#FFC10740"
p_qc$layers[[1]]$geom_params$alpha  <- 0.4
p_qc$layers[[2]]$aes_params$linewidth <- 1.0

t0_qc       <- as.Date(min(p_qc$data$ds))
xint35_qc   <- min(p_qc$data$ds) + lubridate::days(34)

p_qc +
  scale_color_manual(
    values = c(TimeGPT = "red"),
    labels = c(TimeGPT = "Prediction (95% IC in gray area)")
  ) +
  scale_x_datetime(
    name   = "Days of follow-up (35 before, 7 after)",
    limits = c(min(p_qc$data$ds), max(p_qc$data$ds)),
    breaks = seq(min(p_qc$data$ds), max(p_qc$data$ds), by = "7 days"),
    labels = function(x) as.integer(as.Date(x) - t0_qc + 1)
  ) +
  theme_classic(base_size = 13) +
  theme(
    strip.text       = element_blank(),
    strip.background = element_blank(),
    legend.position  = "bottom"
  ) +
  geom_vline(
    xintercept = xint35_qc,
    linetype   = "dashed",
    color      = "black",
    linewidth  = 1
  ) +
  labs(y = "High-speed collisions per 1 M vehicles")
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Prediction on the treated (Québec control only)

Prediction on the treated (Québec control only)

Code
## ── Y_lag10 en un solo pipe ───────────────────────────────────────────────────
collisions_weather_corr_rect_quebec_synth_tr$Y_lag10 <- 
  collisions_weather_corr_rect_quebec_synth_tr |>
    arrange(tr_contr, year.x, yday_corr) |>
    mutate(
      Y_lag10 = dplyr::lag(rate_veh, 10, order_by = yday_corr),
      Y_lag10 = ifelse(is.na(Y_lag10),
                       first(na.omit(Y_lag10)),
                       Y_lag10)
    ) |>
    ungroup() |>
    pull(Y_lag10)

## ── Conjuntos 33 días pre-intervención y 7 post ───────────────────────────────
df33_qc <- collisions_weather_corr_rect_quebec_synth_tr %>% 
  arrange(tr_contr, year.x, yday_corr) %>% 
  filter(yday_corr < 34) %>% 
  transmute(
    unique_id = "treated",
    ds        = as.POSIXct(as.Date("2019-05-19") + (yday_corr - 1)),
    y         = rate_veh,
    mean_median_lag_2_prec_median_imp,
    mean_median_total_precip_median_lin,
    mean_min_temp_mean_lin,
    mean_max_temp_mean_lin,
    Y_lag10
  )

df33_post_qc <- collisions_weather_corr_rect_quebec_synth_tr %>% 
  arrange(tr_contr, year.x, yday_corr) %>% 
  filter(yday_corr >= 34, yday_corr <= max_time) %>% 
  transmute(
    unique_id = "treated",
    ds        = as.POSIXct(as.Date("2019-05-19") + (yday_corr - 1)),
    y         = rate_veh,
    mean_median_lag_2_prec_median_imp,
    mean_median_total_precip_median_lin,
    mean_min_temp_mean_lin,
    mean_max_temp_mean_lin,
    Y_lag10
  )

## ── TimeGPT ───────────────────────────────────────────────────────────────────
fc33_qc <- nixtla_client_forecast(
  rbind(mutate(df33_qc[1:2,], ds = ds - 3), df33_qc),
  h     = 7,
  level = c(95),
  freq  = "D",
  finetune_steps = 10,
  finetune_loss  = "rmse"
)
Using historical exogenous features: mean_median_lag_2_prec_median_imp, mean_median_total_precip_median_lin, mean_min_temp_mean_lin, mean_max_temp_mean_lin, Y_lag10
Code
p33_qc <- nixtla_client_plot(
  rbind(mutate(df33_qc[1:2,], ds = ds - 3),
        df33_qc,
        df33_post_qc),
  fc33_qc,
  max_insample_length = 100
)

## ── Estética idéntica a tu original ───────────────────────────────────────────
p33_qc$layers[[1]]$geom_params$fill   <- "#FFC10740"
p33_qc$layers[[1]]$geom_params$alpha  <- 0.4
p33_qc$layers[[2]]$aes_params$linewidth <- 1.0

t0_33_qc  <- as.Date(min(p33_qc$data$ds))
xint33_qc <- min(p33_qc$data$ds) + lubridate::days(32)

p33_qc +
  scale_color_manual(
    values = c(TimeGPT = "red"),
    labels = c(TimeGPT = "Prediction (95% IC in gray area)")
  ) +
  scale_x_datetime(
    name   = "Days of follow-up (33 before, 7 after)",
    limits = c(min(p33_qc$data$ds), max(p33_qc$data$ds)),
    breaks = seq(min(p33_qc$data$ds), max(p33_qc$data$ds), by = "7 days"),
    labels = function(x) as.integer(as.Date(x) - t0_33_qc + 1)
  ) +
  theme_classic(base_size = 13) +
  theme(
    strip.text       = element_blank(),
    strip.background = element_blank(),
    legend.position  = "bottom"
  ) +
  geom_vline(
    xintercept = xint33_qc,
    linetype   = "dashed",
    color      = "black",
    linewidth  = 1
  ) +
  labs(y = "High-speed collisions per 1 M vehicles")
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Prediction on the treated (-2 days before race, 7 days after — Québec)

Prediction on the treated (-2 days before race, 7 days after — Québec)

Causal impact

We used the CausalImpact package, useful for time series data, inferring a counterfactual post race period based on pre-intervention components1.

Code
# Model 2
ss2d <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ss2d <- AddLocalLevel(ss2d,c(as.numeric(df35$y), rep(NA,7))
                      ) #
# Add weekly seasonal
ss2d <- AddSeasonal(ss2d,c(as.numeric(df35$y), rep(NA,7)), 
                    nseasons=7, season.duration = 1) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
ss2d <- AddAutoAr(ss2d, y = c(as.numeric(df35$y), rep(NA,7)), lags = 10) #NO PUEDO AREGAR AR1 CON POISSON
# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).

# Add external regressors explicitly (e.g., temperatures, precipitation)
covariates <- cbind(
 # y = as.numeric(df35$y),
  mean_min_temp_mean_lin = df35$mean_min_temp_mean_lin,
  mean_max_temp_mean_lin = df35$mean_max_temp_mean_lin,
  mean_median_lag_2_prec_median_imp = df35$mean_median_lag_2_prec_median_imp,
  mean_median_total_precip_median_lin = df35$mean_median_total_precip_median_lin,
  Y_lag10 = df35$Y_lag10
)
# Ensure names explicitly
colnames(covariates) <- c(#"y",
                          "mean_min_temp_mean_lin",
                          "mean_max_temp_mean_lin",
                          "mean_median_lag_2_prec_median_imp",
                          "mean_median_total_precip_median_lin",
                          "Y_lag10")
covariates<- as.data.frame(covariates)

covariates_future <- as.matrix(df35_post[, c("mean_min_temp_mean_lin",
                                             "mean_max_temp_mean_lin",
                                             "mean_median_lag_2_prec_median_imp",
                                             "mean_median_total_precip_median_lin",
                                             "Y_lag10")])

covariates_full <- rbind(covariates, covariates_future)

# Include regression terms
ss2d <- AddDynamicRegression(ss2d, covariates_full)

model2d_ratio <- bsts(c(as.numeric(df35$y), rep(NA,7)),
               state.specification = ss2d, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.
               family ="student", #A Bayesian Analysis of Time-Series Event Count Data
               niter = 3e4, 
               #burn = 2000, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
               seed= 2125)

## Since we have predictor variables in the model, we need to explicitly
# make their coefficients time-varying using AddDynamicRegression().

impact3d_ratio_resp <- CausalImpact(bsts.model = model2d_ratio,model.args = list(prior.level.sd=.1, dynamic.regression=T),
              post.period.response = df35_post$y)

plot(impact3d_ratio_resp, "original") 
Causal impact (Main)

Causal impact (Main)

Code
summary(impact3d_ratio_resp)
plot(impact3d_ratio_resp)
Causal impact (Main)

Causal impact (Main)

=-=-=-=-= Iteration 0 Tue Jul  8 00:27:32 2025
 =-=-=-=-=
=-=-=-=-= Iteration 3000 Tue Jul  8 00:27:35 2025
 =-=-=-=-=
=-=-=-=-= Iteration 6000 Tue Jul  8 00:27:38 2025
 =-=-=-=-=
=-=-=-=-= Iteration 9000 Tue Jul  8 00:27:41 2025
 =-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Jul  8 00:27:44 2025
 =-=-=-=-=
=-=-=-=-= Iteration 15000 Tue Jul  8 00:27:47 2025
 =-=-=-=-=
=-=-=-=-= Iteration 18000 Tue Jul  8 00:27:50 2025
 =-=-=-=-=
=-=-=-=-= Iteration 21000 Tue Jul  8 00:27:53 2025
 =-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Jul  8 00:27:56 2025
 =-=-=-=-=
=-=-=-=-= Iteration 27000 Tue Jul  8 00:27:59 2025
 =-=-=-=-=
Posterior inference {CausalImpact}

                         Average        Cumulative  
Actual                   1.5            10.5        
Prediction (s.d.)        1.6 (0.49)     11.5 (3.41) 
95% CI                   [0.7, 2.6]     [4.9, 18.3] 
                                                    
Absolute effect (s.d.)   -0.14 (0.49)   -0.97 (3.41)
95% CI                   [-1.1, 0.81]   [-7.8, 5.65]
                                                    
Relative effect (s.d.)   -1.3% (766%)   -1.3% (766%)
95% CI                   [-43%, 113%]   [-43%, 113%]

Posterior tail-area probability p:   0.38543
Posterior prob. of a causal effect:  61%

For more details, type: summary(impact, "report")

We included other terms (local linear trend).

Code
library(bsts)

y35 <- c(as.numeric(df35$y), rep(NA, 7)) # respuesta escalada por 100 millones

ss2e <- list()
ss2e <- AddLocalLinearTrend(ss2e, y35)
ss2e <- AddSeasonal(ss2e, y35, nseasons = 7, season.duration = 1)
ss2e <- AddAutoAr(ss2e, y35, lags = 1)


ss2e <- AddDynamicRegression(ss2e, covariates_full)

set.seed(2125)
model35 <- bsts(y35, state.specification = ss2e, seed=2125,
              niter = 3e4, family = "student")

#plot(model, "components", burn = 2000)

impact3d_ratio_resp2 <- CausalImpact(bsts.model = model35,model.args = list(prior.level.sd=.1, dynamic.regression=T),
              post.period.response = df35_post$y)
#plot(impact3d, "original") 

summary(impact3d_ratio_resp2)

plot(impact3d_ratio_resp2)
Causal impact (Alternative, w/ linear local trend)

Causal impact (Alternative, w/ linear local trend)

=-=-=-=-= Iteration 0 Tue Jul  8 00:28:10 2025
 =-=-=-=-=
=-=-=-=-= Iteration 3000 Tue Jul  8 00:28:13 2025
 =-=-=-=-=
=-=-=-=-= Iteration 6000 Tue Jul  8 00:28:15 2025
 =-=-=-=-=
=-=-=-=-= Iteration 9000 Tue Jul  8 00:28:18 2025
 =-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Jul  8 00:28:21 2025
 =-=-=-=-=
=-=-=-=-= Iteration 15000 Tue Jul  8 00:28:23 2025
 =-=-=-=-=
=-=-=-=-= Iteration 18000 Tue Jul  8 00:28:26 2025
 =-=-=-=-=
=-=-=-=-= Iteration 21000 Tue Jul  8 00:28:28 2025
 =-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Jul  8 00:28:30 2025
 =-=-=-=-=
=-=-=-=-= Iteration 27000 Tue Jul  8 00:28:33 2025
 =-=-=-=-=
Posterior inference {CausalImpact}

                         Average        Cumulative  
Actual                   1.5            10.5        
Prediction (s.d.)        1.7 (0.53)     11.8 (3.70) 
95% CI                   [0.67, 2.8]    [4.69, 19.5]
                                                    
Absolute effect (s.d.)   -0.19 (0.53)   -1.30 (3.70)
95% CI                   [-1.3, 0.83]   [-8.9, 5.83]
                                                    
Relative effect (s.d.)   12% (2224%)    12% (2224%) 
95% CI                   [-46%, 119%]   [-46%, 119%]

Posterior tail-area probability p:   0.35602
Posterior prob. of a causal effect:  64%

For more details, type: summary(impact, "report")
-2 days pre
Code
# Model 2
ss2a <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ss2a <- AddLocalLevel(ss2a,c(as.numeric(df33$y), rep(NA,7))
                      ) #
# Add weekly seasonal
ss2a <- AddSeasonal(ss2a,c(as.numeric(df33$y), rep(NA,7)), 
                    nseasons=7, season.duration = 1) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
ss2a <- AddAutoAr(ss2a, y = c(as.numeric(df33$y), rep(NA,7)), lags = 10) #NO PUEDO AREGAR AR1 CON POISSON
# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).

# Add external regressors explicitly (e.g., temperatures, precipitation)
covariates33 <- cbind(
 # y = as.numeric(df35$y),
  mean_min_temp_mean_lin = df33$mean_min_temp_mean_lin,
  mean_max_temp_mean_lin = df33$mean_max_temp_mean_lin,
  mean_median_lag_2_prec_median_imp = df33$mean_median_lag_2_prec_median_imp,
  mean_median_total_precip_median_lin= df33$mean_median_total_precip_median_lin,
  Y_lag10 = df33$Y_lag10
)
# Ensure names explicitly
colnames(covariates33) <- c(#"y",
                          "mean_min_temp_mean_lin",
                          "mean_max_temp_mean_lin",
                          "mean_median_lag_2_prec_median_imp",
                          "mean_median_total_precip_median_lin",
                          "Y_lag10")
covariates33<- as.data.frame(covariates33)

covariates_future33 <- as.matrix(df33_post[1:7, c("mean_min_temp_mean_lin",
                                             "mean_max_temp_mean_lin",
                                             "mean_median_lag_2_prec_median_imp",
                                             "mean_median_total_precip_median_lin",
                                             "Y_lag10")])

covariates_full33 <- rbind(covariates33, covariates_future33)

# Include regression terms
ss2a <- AddDynamicRegression(ss2a, covariates_full33)

model2a_ratio <- bsts(c(as.numeric(df33$y), rep(NA,7)),
               state.specification = ss2a, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.
               family ="student", #A Bayesian Analysis of Time-Series Event Count Data
               niter = 3e4, 
               #burn = 2000, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
               seed= 2125)

## Since we have predictor variables in the model, we need to explicitly
# make their coefficients time-varying using AddDynamicRegression().

impact3a_ratio_resp <- CausalImpact(bsts.model = model2a_ratio,model.args = list(prior.level.sd=.1, dynamic.regression=T),
              post.period.response = df33_post$y[1:7])

plot(impact3a_ratio_resp, "original") 
Causal impact (Main, -2 days race)

Causal impact (Main, -2 days race)

Code
summary(impact3a_ratio_resp)
plot(impact3a_ratio_resp)
Causal impact (Main, -2 days race)

Causal impact (Main, -2 days race)

=-=-=-=-= Iteration 0 Tue Jul  8 00:28:40 2025
 =-=-=-=-=
=-=-=-=-= Iteration 3000 Tue Jul  8 00:28:43 2025
 =-=-=-=-=
=-=-=-=-= Iteration 6000 Tue Jul  8 00:28:46 2025
 =-=-=-=-=
=-=-=-=-= Iteration 9000 Tue Jul  8 00:28:49 2025
 =-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Jul  8 00:28:52 2025
 =-=-=-=-=
=-=-=-=-= Iteration 15000 Tue Jul  8 00:28:55 2025
 =-=-=-=-=
=-=-=-=-= Iteration 18000 Tue Jul  8 00:28:58 2025
 =-=-=-=-=
=-=-=-=-= Iteration 21000 Tue Jul  8 00:29:01 2025
 =-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Jul  8 00:29:04 2025
 =-=-=-=-=
=-=-=-=-= Iteration 27000 Tue Jul  8 00:29:07 2025
 =-=-=-=-=
Posterior inference {CausalImpact}

                         Average         Cumulative   
Actual                   1.5             10.4         
Prediction (s.d.)        1.6 (0.43)      11.0 (2.98)  
95% CI                   [0.74, 2.4]     [5.19, 16.9] 
                                                      
Absolute effect (s.d.)   -0.086 (0.43)   -0.599 (2.98)
95% CI                   [-0.94, 0.74]   [-6.55, 5.17]
                                                      
Relative effect (s.d.)   2.2% (258%)     2.2% (258%)  
95% CI                   [-39%, 99%]     [-39%, 99%]  

Posterior tail-area probability p:   0.4221
Posterior prob. of a causal effect:  58%

For more details, type: summary(impact, "report")
Code
library(bsts)

y33b <- c(as.numeric(df33$y), rep(NA, 7)) # respuesta escalada por 100 millones

ss2b <- list()
ss2b <- AddLocalLinearTrend(ss2b, y33b)
ss2b <- AddSeasonal(ss2b, y33b, nseasons = 7, season.duration = 1)
ss2b <- AddAutoAr(ss2b, y33b, lags = 1)


ss2b <- AddDynamicRegression(ss2b, covariates_full33)

set.seed(2125)
model33b <- bsts(y33b, state.specification = ss2b, seed=2125,
              niter = 3e4, family = "student")

#plot(model, "components", burn = 2000)

impact3d_ratio_resp2 <- CausalImpact(bsts.model = model33b,model.args = list(prior.level.sd=.1, dynamic.regression=T),
              post.period.response = df33_post$y[1:7])
#plot(impact3d, "original") 

summary(impact3d_ratio_resp2)

plot(impact3d_ratio_resp2)
Causal impact (Local linear trend, -2 days race)

Causal impact (Local linear trend, -2 days race)

=-=-=-=-= Iteration 0 Tue Jul  8 00:29:17 2025
 =-=-=-=-=
=-=-=-=-= Iteration 3000 Tue Jul  8 00:29:20 2025
 =-=-=-=-=
=-=-=-=-= Iteration 6000 Tue Jul  8 00:29:22 2025
 =-=-=-=-=
=-=-=-=-= Iteration 9000 Tue Jul  8 00:29:25 2025
 =-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Jul  8 00:29:27 2025
 =-=-=-=-=
=-=-=-=-= Iteration 15000 Tue Jul  8 00:29:30 2025
 =-=-=-=-=
=-=-=-=-= Iteration 18000 Tue Jul  8 00:29:32 2025
 =-=-=-=-=
=-=-=-=-= Iteration 21000 Tue Jul  8 00:29:34 2025
 =-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Jul  8 00:29:37 2025
 =-=-=-=-=
=-=-=-=-= Iteration 27000 Tue Jul  8 00:29:40 2025
 =-=-=-=-=
Posterior inference {CausalImpact}

                         Average        Cumulative  
Actual                   1.5            10.4        
Prediction (s.d.)        1.6 (0.47)     11.4 (3.27) 
95% CI                   [0.74, 2.6]    [5.17, 18.1]
                                                    
Absolute effect (s.d.)   -0.15 (0.47)   -1.06 (3.27)
95% CI                   [-1.1, 0.74]   [-7.7, 5.20]
                                                    
Relative effect (s.d.)   -6.3% (674%)   -6.3% (674%)
95% CI                   [-43%, 99%]    [-43%, 99%] 

Posterior tail-area probability p:   0.37243
Posterior prob. of a causal effect:  63%

For more details, type: summary(impact, "report")
-2 days pre, +3
Code
# Model 2
ss2ab <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ss2a <- AddLocalLevel(ss2ab,c(as.numeric(df33$y), rep(NA,3))
                      ) #
# Add weekly seasonal
ss2ab <- AddSeasonal(ss2ab,c(as.numeric(df33$y), rep(NA,3)), 
                    nseasons=7, season.duration = 1) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
ss2ab <- AddAutoAr(ss2ab, y = c(as.numeric(df33$y), rep(NA,3)), lags = 10) #NO PUEDO AREGAR AR1 CON POISSON
# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).

covariates_future33b <- as.matrix(df33_post[1:3, c("mean_min_temp_mean_lin",
                                             "mean_max_temp_mean_lin",
                                             "mean_median_lag_2_prec_median_imp",
                                             "mean_median_total_precip_median_lin",
                                             "Y_lag10")])

covariates_full33b <- rbind(covariates33, covariates_future33b)

# Include regression terms
ss2ab <- AddDynamicRegression(ss2ab, covariates_full33b)

model2ab_ratio <- bsts(c(as.numeric(df33$y), rep(NA,3)),
               state.specification = ss2ab, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.
               family ="student", #A Bayesian Analysis of Time-Series Event Count Data
               niter = 3e4, 
               #burn = 2000, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
               seed= 2125)

## Since we have predictor variables in the model, we need to explicitly
# make their coefficients time-varying using AddDynamicRegression().

impact3ab_ratio_resp <- CausalImpact(bsts.model = model2ab_ratio,model.args = list(prior.level.sd=.1, dynamic.regression=T),
              post.period.response = df33_post$y[1:3])

plot(impact3ab_ratio_resp, "original") 
Causal impact (Main, -2 days race, +3 days post)

Causal impact (Main, -2 days race, +3 days post)

Code
summary(impact3ab_ratio_resp)
plot(impact3ab_ratio_resp)
Causal impact (Main, -2 days race, +3 days post)

Causal impact (Main, -2 days race, +3 days post)

=-=-=-=-= Iteration 0 Tue Jul  8 00:29:47 2025
 =-=-=-=-=
=-=-=-=-= Iteration 3000 Tue Jul  8 00:29:49 2025
 =-=-=-=-=
=-=-=-=-= Iteration 6000 Tue Jul  8 00:29:52 2025
 =-=-=-=-=
=-=-=-=-= Iteration 9000 Tue Jul  8 00:29:54 2025
 =-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Jul  8 00:29:57 2025
 =-=-=-=-=
=-=-=-=-= Iteration 15000 Tue Jul  8 00:29:59 2025
 =-=-=-=-=
=-=-=-=-= Iteration 18000 Tue Jul  8 00:30:01 2025
 =-=-=-=-=
=-=-=-=-= Iteration 21000 Tue Jul  8 00:30:04 2025
 =-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Jul  8 00:30:06 2025
 =-=-=-=-=
=-=-=-=-= Iteration 27000 Tue Jul  8 00:30:09 2025
 =-=-=-=-=
Posterior inference {CausalImpact}

                         Average         Cumulative   
Actual                   1.6             4.8          
Prediction (s.d.)        1.7 (0.35)      5.1 (1.04)   
95% CI                   [1.0, 2.4]      [3.1, 7.2]   
                                                      
Absolute effect (s.d.)   -0.11 (0.35)    -0.34 (1.04) 
95% CI                   [-0.81, 0.56]   [-2.42, 1.68]
                                                      
Relative effect (s.d.)   -2.1% (26%)     -2.1% (26%)  
95% CI                   [-34%, 54%]     [-34%, 54%]  

Posterior tail-area probability p:   0.37121
Posterior prob. of a causal effect:  63%

For more details, type: summary(impact, "report")
Code
library(bsts)

y33b <- c(as.numeric(df33$y), rep(NA, 3)) # respuesta escalada por 100 millones

ss2b <- list()
ss2b <- AddLocalLinearTrend(ss2b, y33b)
ss2b <- AddSeasonal(ss2b, y33b, nseasons = 7, season.duration = 1)
ss2b <- AddAutoAr(ss2b, y33b, lags = 1)


ss2b <- AddDynamicRegression(ss2b, covariates_full33b)

set.seed(2125)
model33b <- bsts(y33b, state.specification = ss2b, seed=2125,
              niter = 3e4, family = "student")

#plot(model, "components", burn = 2000)

impact3db_ratio_resp2 <- CausalImpact(bsts.model = model33b, model.args = list(prior.level.sd=.1, dynamic.regression=T),
              post.period.response = df33_post$y[1:3])
#plot(impact3d, "original") 

summary(impact3db_ratio_resp2)

plot(impact3db_ratio_resp2)
Causal impact (Local linear trend, -2 days race, 3)

Causal impact (Local linear trend, -2 days race, 3)

=-=-=-=-= Iteration 0 Tue Jul  8 00:30:19 2025
 =-=-=-=-=
=-=-=-=-= Iteration 3000 Tue Jul  8 00:30:21 2025
 =-=-=-=-=
=-=-=-=-= Iteration 6000 Tue Jul  8 00:30:23 2025
 =-=-=-=-=
=-=-=-=-= Iteration 9000 Tue Jul  8 00:30:26 2025
 =-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Jul  8 00:30:28 2025
 =-=-=-=-=
=-=-=-=-= Iteration 15000 Tue Jul  8 00:30:30 2025
 =-=-=-=-=
=-=-=-=-= Iteration 18000 Tue Jul  8 00:30:32 2025
 =-=-=-=-=
=-=-=-=-= Iteration 21000 Tue Jul  8 00:30:34 2025
 =-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Jul  8 00:30:37 2025
 =-=-=-=-=
=-=-=-=-= Iteration 27000 Tue Jul  8 00:30:39 2025
 =-=-=-=-=
Posterior inference {CausalImpact}

                         Average         Cumulative   
Actual                   1.6             4.8          
Prediction (s.d.)        1.6 (0.36)      4.9 (1.07)   
95% CI                   [0.95, 2.4]     [2.84, 7.1]  
                                                      
Absolute effect (s.d.)   -0.043 (0.36)   -0.130 (1.07)
95% CI                   [-0.76, 0.65]   [-2.28, 1.95]
                                                      
Relative effect (s.d.)   3% (30%)        3% (30%)     
95% CI                   [-32%, 69%]     [-32%, 69%]  

Posterior tail-area probability p:   0.44954
Posterior prob. of a causal effect:  55%

For more details, type: summary(impact, "report")

Only Quebec as controls

Code
## ── Estado: nivel local + estacionalidad + AR(10) ────────────────────────────
ss2d_qc <- list()
ss2d_qc <- AddLocalLevel(
  ss2d_qc, c(as.numeric(df35_qc$y), rep(NA, 7))
)
ss2d_qc <- AddSeasonal(
  ss2d_qc, c(as.numeric(df35_qc$y), rep(NA, 7)),
  nseasons = 7, season.duration = 1
)
ss2d_qc <- AddAutoAr(
  ss2d_qc, y = c(as.numeric(df35_qc$y), rep(NA, 7)),
  lags = 10
)

## ── Covariables como data.frame ──────────────────────────────────────────────
covariates_qc <- df35_qc |> dplyr::select(
  mean_min_temp_mean_lin,
  mean_max_temp_mean_lin,
  mean_median_lag_2_prec_median_imp,
  mean_median_total_precip_median_lin,
  Y_lag10
)

covariates_future_qc <- df35_post_qc |> dplyr::select(
  mean_min_temp_mean_lin,
  mean_max_temp_mean_lin,
  mean_median_lag_2_prec_median_imp,
  mean_median_total_precip_median_lin,
  Y_lag10
)

covariates_full_qc <- dplyr::bind_rows(covariates_qc, covariates_future_qc)

## ── Regresión dinámica ───────────────────────────────────────────────────────
ss2d_qc <- AddDynamicRegression(ss2d_qc, covariates_full_qc)

## ── Ajuste BSTS y CausalImpact ───────────────────────────────────────────────
model2d_qc <- bsts(
  c(as.numeric(df35_qc$y), rep(NA, 7)),
  state.specification = ss2d_qc,
  family = "student",
  niter  = 3e4,
  seed   = 2125
)

impact_qc <- CausalImpact(
  bsts.model           = model2d_qc,
  model.args           = list(prior.level.sd = .1, dynamic.regression = TRUE),
  post.period.response = df35_post_qc$y
)

plot(impact_qc, "original")
Causal impact (Québec only, 0 days pre)

Causal impact (Québec only, 0 days pre)

Code
summary(impact_qc)
=-=-=-=-= Iteration 0 Tue Jul  8 00:30:45 2025
 =-=-=-=-=
=-=-=-=-= Iteration 3000 Tue Jul  8 00:30:48 2025
 =-=-=-=-=
=-=-=-=-= Iteration 6000 Tue Jul  8 00:30:51 2025
 =-=-=-=-=
=-=-=-=-= Iteration 9000 Tue Jul  8 00:30:54 2025
 =-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Jul  8 00:30:57 2025
 =-=-=-=-=
=-=-=-=-= Iteration 15000 Tue Jul  8 00:31:00 2025
 =-=-=-=-=
=-=-=-=-= Iteration 18000 Tue Jul  8 00:31:04 2025
 =-=-=-=-=
=-=-=-=-= Iteration 21000 Tue Jul  8 00:31:07 2025
 =-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Jul  8 00:31:10 2025
 =-=-=-=-=
=-=-=-=-= Iteration 27000 Tue Jul  8 00:31:13 2025
 =-=-=-=-=
Posterior inference {CausalImpact}

                         Average         Cumulative   
Actual                   1.4             10.0         
Prediction (s.d.)        1.5 (0.46)      10.6 (3.23)  
95% CI                   [0.62, 2.4]     [4.35, 17.1] 
                                                      
Absolute effect (s.d.)   -0.079 (0.46)   -0.554 (3.23)
95% CI                   [-1.0, 0.81]    [-7.1, 5.70] 
                                                      
Relative effect (s.d.)   6.2% (216%)     6.2% (216%)  
95% CI                   [-42%, 126%]    [-42%, 126%] 

Posterior tail-area probability p:   0.43895
Posterior prob. of a causal effect:  56%

For more details, type: summary(impact, "report")
Code
## ── Vector de respuesta ───────────────────────────────────────────────────────
y35_qc <- c(as.numeric(df35_qc$y), rep(NA, 7))

## ── Covariables como data.frame ──────────────────────────────────────────────
covariates_qc <- df35_qc |> dplyr::select(
  mean_min_temp_mean_lin,
  mean_max_temp_mean_lin,
  mean_median_lag_2_prec_median_imp,
  mean_median_total_precip_median_lin,
  Y_lag10
)

covariates_future_qc <- df35_post_qc |> dplyr::select(
  mean_min_temp_mean_lin,
  mean_max_temp_mean_lin,
  mean_median_lag_2_prec_median_imp,
  mean_median_total_precip_median_lin,
  Y_lag10
)

covariates_full_qc <- dplyr::bind_rows(covariates_qc, covariates_future_qc)

## ── Especificación del estado: tendencia local lineal ────────────────────────
ss2e_qc <- list()
ss2e_qc <- AddLocalLinearTrend(ss2e_qc, y35_qc)
ss2e_qc <- AddSeasonal(ss2e_qc, y35_qc, nseasons = 7, season.duration = 1)
ss2e_qc <- AddAutoAr(ss2e_qc, y35_qc, lags = 1)
ss2e_qc <- AddDynamicRegression(ss2e_qc, covariates_full_qc)

## ── Ajuste BSTS y CausalImpact ───────────────────────────────────────────────
model35_qc <- bsts(
  y35_qc,
  state.specification = ss2e_qc,
  family = "student",
  niter  = 3e4,
  seed   = 2125
)

impact2_qc <- CausalImpact(
  bsts.model           = model35_qc,
  model.args           = list(prior.level.sd = .1, dynamic.regression = TRUE),
  post.period.response = df35_post_qc$y
)

summary(impact2_qc)
plot(impact2_qc)
Causal impact (Québec only, 0 days pre, local linear trend)

Causal impact (Québec only, 0 days pre, local linear trend)

=-=-=-=-= Iteration 0 Tue Jul  8 00:31:20 2025
 =-=-=-=-=
=-=-=-=-= Iteration 3000 Tue Jul  8 00:31:22 2025
 =-=-=-=-=
=-=-=-=-= Iteration 6000 Tue Jul  8 00:31:25 2025
 =-=-=-=-=
=-=-=-=-= Iteration 9000 Tue Jul  8 00:31:27 2025
 =-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Jul  8 00:31:30 2025
 =-=-=-=-=
=-=-=-=-= Iteration 15000 Tue Jul  8 00:31:32 2025
 =-=-=-=-=
=-=-=-=-= Iteration 18000 Tue Jul  8 00:31:35 2025
 =-=-=-=-=
=-=-=-=-= Iteration 21000 Tue Jul  8 00:31:37 2025
 =-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Jul  8 00:31:40 2025
 =-=-=-=-=
=-=-=-=-= Iteration 27000 Tue Jul  8 00:31:42 2025
 =-=-=-=-=
Posterior inference {CausalImpact}

                         Average          Cumulative    
Actual                   1.4              10.0          
Prediction (s.d.)        1.5 (0.5)        10.8 (3.5)    
95% CI                   [0.56, 2.6]      [3.94, 17.9]  
                                                        
Absolute effect (s.d.)   -0.11 (0.5)      -0.76 (3.5)   
95% CI                   [-1.1, 0.87]     [-7.9, 6.10]  
                                                        
Relative effect (s.d.)   -0.74% (1070%)   -0.74% (1070%)
95% CI                   [-44%, 147%]     [-44%, 147%]  

Posterior tail-area probability p:   0.41321
Posterior prob. of a causal effect:  59%

For more details, type: summary(impact, "report")
-2 days pre
Code
## ── Vector de respuesta ───────────────────────────────────────────────────────
y33_qc <- c(as.numeric(df33_qc$y), rep(NA, 7))  # 7 post-días + 2 de ajuste

## ── Covariables (data.frame) ─────────────────────────────────────────────────
covariates33_qc <- df33_qc |> dplyr::select(
  mean_min_temp_mean_lin,
  mean_max_temp_mean_lin,
  mean_median_lag_2_prec_median_imp,
  mean_median_total_precip_median_lin,
  Y_lag10
)

covariates_future33_qc <- df33_post_qc |> dplyr::select(
  mean_min_temp_mean_lin,
  mean_max_temp_mean_lin,
  mean_median_lag_2_prec_median_imp,
  mean_median_total_precip_median_lin,
  Y_lag10
)

covariates_full33_qc <- dplyr::bind_rows(covariates33_qc, covariates_future33_qc[1:7,])

## ── Especificación del estado: nivel local ───────────────────────────────────
ss2a_qc <- list()
ss2a_qc <- AddLocalLevel(ss2a_qc, y33_qc)
ss2a_qc <- AddSeasonal(ss2a_qc, y33_qc, nseasons = 7, season.duration = 1)
ss2a_qc <- AddAutoAr   (ss2a_qc, y33_qc, lags = 10)
ss2a_qc <- AddDynamicRegression(ss2a_qc, covariates_full33_qc)

## ── Ajuste BSTS y CausalImpact ───────────────────────────────────────────────
model2a_qc <- bsts(
  y33_qc,
  state.specification = ss2a_qc,
  family = "student",
  niter  = 3e4,
  seed   = 2125
)

impact33_qc <- CausalImpact(
  bsts.model           = model2a_qc,
  model.args           = list(prior.level.sd = .1, dynamic.regression = TRUE),
  post.period.response = df33_post_qc$y[1:7]
)

plot(impact33_qc, "original")
Causal impact (Québec only -2 days pre, model w/ local level)

Causal impact (Québec only -2 days pre, model w/ local level)

Code
summary(impact33_qc)
=-=-=-=-= Iteration 0 Tue Jul  8 00:31:49 2025
 =-=-=-=-=
=-=-=-=-= Iteration 3000 Tue Jul  8 00:31:52 2025
 =-=-=-=-=
=-=-=-=-= Iteration 6000 Tue Jul  8 00:31:55 2025
 =-=-=-=-=
=-=-=-=-= Iteration 9000 Tue Jul  8 00:31:58 2025
 =-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Jul  8 00:32:01 2025
 =-=-=-=-=
=-=-=-=-= Iteration 15000 Tue Jul  8 00:32:04 2025
 =-=-=-=-=
=-=-=-=-= Iteration 18000 Tue Jul  8 00:32:07 2025
 =-=-=-=-=
=-=-=-=-= Iteration 21000 Tue Jul  8 00:32:10 2025
 =-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Jul  8 00:32:13 2025
 =-=-=-=-=
=-=-=-=-= Iteration 27000 Tue Jul  8 00:32:16 2025
 =-=-=-=-=
Posterior inference {CausalImpact}

                         Average         Cumulative   
Actual                   1.4             9.8          
Prediction (s.d.)        1.5 (0.4)       10.2 (2.8)   
95% CI                   [0.68, 2.3]     [4.76, 15.8] 
                                                      
Absolute effect (s.d.)   -0.066 (0.4)    -0.460 (2.8) 
95% CI                   [-0.86, 0.72]   [-6.03, 5.03]
                                                      
Relative effect (s.d.)   -3.9% (1355%)   -3.9% (1355%)
95% CI                   [-38%, 104%]    [-38%, 104%] 

Posterior tail-area probability p:   0.43406
Posterior prob. of a causal effect:  57%

For more details, type: summary(impact, "report")
Code
## ── Vector de respuesta ───────────────────────────────────────────────────────
y33_lin_qc <- c(as.numeric(df33_qc$y), rep(NA, 7))

## ── Covariables (data.frame) ─────────────────────────────────────────────────
covariates_full33_qc <- dplyr::bind_rows(covariates33_qc, covariates_future33_qc)

## ── Estado: tendencia local lineal ───────────────────────────────────────────
ss2b_qc <- list()
ss2b_qc <- AddLocalLinearTrend(ss2b_qc, y33_lin_qc)
ss2b_qc <- AddSeasonal        (ss2b_qc, y33_lin_qc, nseasons = 7, season.duration = 1)
ss2b_qc <- AddAutoAr          (ss2b_qc, y33_lin_qc, lags = 1)
ss2b_qc <- AddDynamicRegression(ss2b_qc, covariates_full33_qc)

## ── Ajuste BSTS y CausalImpact ───────────────────────────────────────────────
model33b_qc <- bsts(
  y33_lin_qc,
  state.specification = ss2b_qc,
  family = "student",
  niter  = 3e4,
  seed   = 2125
)

impact2_33_qc <- CausalImpact(
  bsts.model           = model33b_qc,
  model.args           = list(prior.level.sd = .1, dynamic.regression = TRUE),
  post.period.response = df33_post_qc$y[1:7]
)

summary(impact2_33_qc)
plot(impact2_33_qc)
Causal impact (Québec only,-2 days pre, local linear trend)

Causal impact (Québec only,-2 days pre, local linear trend)

=-=-=-=-= Iteration 0 Tue Jul  8 00:32:23 2025
 =-=-=-=-=
=-=-=-=-= Iteration 3000 Tue Jul  8 00:32:25 2025
 =-=-=-=-=
=-=-=-=-= Iteration 6000 Tue Jul  8 00:32:28 2025
 =-=-=-=-=
=-=-=-=-= Iteration 9000 Tue Jul  8 00:32:30 2025
 =-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Jul  8 00:32:32 2025
 =-=-=-=-=
=-=-=-=-= Iteration 15000 Tue Jul  8 00:32:35 2025
 =-=-=-=-=
=-=-=-=-= Iteration 18000 Tue Jul  8 00:32:37 2025
 =-=-=-=-=
=-=-=-=-= Iteration 21000 Tue Jul  8 00:32:39 2025
 =-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Jul  8 00:32:42 2025
 =-=-=-=-=
=-=-=-=-= Iteration 27000 Tue Jul  8 00:32:44 2025
 =-=-=-=-=
Posterior inference {CausalImpact}

                         Average         Cumulative   
Actual                   1.4             9.8          
Prediction (s.d.)        1.5 (0.43)      10.7 (3.00)  
95% CI                   [0.68, 2.4]     [4.73, 16.7] 
                                                      
Absolute effect (s.d.)   -0.12 (0.43)    -0.86 (3.00) 
95% CI                   [-0.98, 0.72]   [-6.87, 5.06]
                                                      
Relative effect (s.d.)   3.9% (164%)     3.9% (164%)  
95% CI                   [-41%, 105%]    [-41%, 105%] 

Posterior tail-area probability p:   0.38691
Posterior prob. of a causal effect:  61%

For more details, type: summary(impact, "report")
-2 days pre, +3
Code
## ── Vector de respuesta ───────────────────────────────────────────────────────
y33b_qc <- c(as.numeric(df33_qc$y), rep(NA, 3))  # 7 post-días + 2 de ajuste

covariates_full33b_qc <- dplyr::bind_rows(covariates33_qc, covariates_future33_qc[1:3,])

## ── Especificación del estado: nivel local ───────────────────────────────────
ss2ab_qc <- list()
ss2ab_qc <- AddLocalLevel(ss2ab_qc, y33b_qc)
ss2ab_qc <- AddSeasonal(ss2ab_qc, y33b_qc, nseasons = 7, season.duration = 1)
ss2ab_qc <- AddAutoAr   (ss2ab_qc, y33b_qc, lags = 10)
ss2ab_qc <- AddDynamicRegression(ss2ab_qc, covariates_full33b_qc)

## ── Ajuste BSTS y CausalImpact ───────────────────────────────────────────────
model2ab_qc <- bsts(
  y33_qc,
  state.specification = ss2ab_qc,
  family = "student",
  niter  = 3e4,
  seed   = 2125
)
Error in bsts(y33_qc, state.specification = ss2ab_qc, family = "student", : Caught exception with the following error message: 
incompatible vector in SparseVector dot product: 
dense vector: 0.0359117 3.55378e-06 -4.29923e-06 9.23099e-08 1.82078e-07 1.1973e-07 4.18066e-07 0 0 0 0 3.27532e-06 -6.80385e-06 4.98775e-06 -2.6448e-06 -1.53716e-06 9.6989e-07 -0.00153241 -0.000415958 -0.000115887 -0.000321038 
sparse[0] = 1
sparse[1] = 1
sparse[7] = 1
Code
impact33b_qc <- CausalImpact(
  bsts.model           = model2ab_qc,
  model.args           = list(prior.level.sd = .1, dynamic.regression = TRUE),
  post.period.response = df33_post_qc$y[1:3]
)
Error: objeto 'model2ab_qc' no encontrado
Code
plot(impact33b_qc, "original")
Error: objeto 'impact33b_qc' no encontrado
Code
summary(impact33b_qc)
Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': objeto 'impact33b_qc' no encontrado
Code
## ── Vector de respuesta ───────────────────────────────────────────────────────
y33b_lin_qc <- c(as.numeric(df33_qc$y), rep(NA, 3))

## ── Covariables (data.frame) ─────────────────────────────────────────────────
covariates_full33b_qc <- dplyr::bind_rows(covariates33_qc, covariates_future33b_qc)
Error: objeto 'covariates_future33b_qc' no encontrado
Code
## ── Estado: tendencia local lineal ───────────────────────────────────────────
ss2bb_qc <- list()
ss2bb_qc <- AddLocalLinearTrend(ss2bb_qc, y33b_lin_qc)
ss2bb_qc <- AddSeasonal        (ss2bb_qc, y33b_lin_qc, nseasons = 7, season.duration = 1)
ss2bb_qc <- AddAutoAr          (ss2bb_qc, y33b_lin_qc, lags = 1)
ss2bb_qc <- AddDynamicRegression(ss2bb_qc, covariates_full33b_qc)

## ── Ajuste BSTS y CausalImpact ───────────────────────────────────────────────
model33b_qc <- bsts(
  y33b_lin_qc,
  state.specification = ss2bb_qc,
  family = "student",
  niter  = 3e4,
  seed   = 2125
)

impact2_33b_qc <- CausalImpact(
  bsts.model           = model33bb_qc,
  model.args           = list(prior.level.sd = .1, dynamic.regression = TRUE),
  post.period.response = df33_post_qc$y[1:3]
)
Error: objeto 'model33bb_qc' no encontrado
Code
summary(impact2_33b_qc)
Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': objeto 'impact2_33b_qc' no encontrado
Code
plot(impact2_33b_qc)
Error: objeto 'impact2_33b_qc' no encontrado
=-=-=-=-= Iteration 0 Tue Jul  8 00:32:52 2025
 =-=-=-=-=
=-=-=-=-= Iteration 3000 Tue Jul  8 00:32:54 2025
 =-=-=-=-=
=-=-=-=-= Iteration 6000 Tue Jul  8 00:32:56 2025
 =-=-=-=-=
=-=-=-=-= Iteration 9000 Tue Jul  8 00:32:58 2025
 =-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Jul  8 00:33:00 2025
 =-=-=-=-=
=-=-=-=-= Iteration 15000 Tue Jul  8 00:33:03 2025
 =-=-=-=-=
=-=-=-=-= Iteration 18000 Tue Jul  8 00:33:05 2025
 =-=-=-=-=
=-=-=-=-= Iteration 21000 Tue Jul  8 00:33:07 2025
 =-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Jul  8 00:33:09 2025
 =-=-=-=-=
=-=-=-=-= Iteration 27000 Tue Jul  8 00:33:11 2025
 =-=-=-=-=


Session info

Code
cat(paste0("R library: ", Sys.getenv("R_LIBS_USER")))
cat(paste0("Date: ",withr::with_locale(new = c('LC_TIME' = 'C'), code =Sys.time())))
cat(paste0("Editor context: ", getwd()))
cat("quarto version: "); system("quarto --version") 

quarto::quarto_version()

save.image("_data/step3.RData")
R library: H:/Mi unidad/PERSONAL ANDRES/UCH_salud_publica/pasantia/f1/f1/renv/library/windows/R-4.4/x86_64-w64-mingw32Date: 2025-07-08 00:33:14.740894Editor context: H:/Mi unidad/PERSONAL ANDRES/UCH_salud_publica/pasantia/f1/f1quarto version: [1] 0
[1] '1.6.39'
Code
sesion_info <- devtools::session_info()
Warning in system2("quarto", "-V", stdout = TRUE, env = paste0("TMPDIR=", : el
comando ejecutado '"quarto"
TMPDIR=C:/Users/andre/AppData/Local/Temp/RtmpYrN5xg/file569457bd5dcb -V' tiene
el estatus 1
Code
dplyr::select(
  tibble::as_tibble(sesion_info$packages),
  c(package, loadedversion, source)
) |> 
 knitr::kable(caption = "R packages", format = "html",
      col.names = c("Row number", "Package", "Version"),
    row.names = FALSE,
      align = c("c", "l", "r")) |> 
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12)|> 
  kableExtra::scroll_box(width = "100%", height = "375px")  
R packages
Row number Package Version
abind 1.4-8 RSPM
assertthat 0.2.1 RSPM
backports 1.5.0 RSPM
bayesplot 1.12.0 RSPM
bdsmatrix 1.3-7 RSPM
bit 4.6.0 CRAN (R 4.4.3)
bit64 4.6.0-1 CRAN (R 4.4.3)
Boom 0.9.15 RSPM
BoomSpikeSlab 1.2.6 RSPM
boot 1.3-30 CRAN (R 4.4.1)
bpmn 0.1.0 Github (bergant/bpmn@628d3efaa27544c221b2fd7f1895301a63b70c49)
bridgesampling 1.1-2 RSPM
brms 2.22.0 RSPM
Brobdingnag 1.2-9 RSPM
broom 1.0.8 RSPM
bsts 0.9.10 RSPM
cachem 1.1.0 CRAN (R 4.4.3)
car 3.1-3 RSPM
carData 3.0-5 RSPM
CausalImpact 1.3.0 RSPM
CBPS 0.23 RSPM
checkmate 2.3.2 RSPM
cli 3.6.5 RSPM
coda 0.19-4.1 RSPM
codetools 0.2-20 CRAN (R 4.4.1)
collapse 2.1.2 RSPM
colorspace 2.1-1 RSPM
conquer 1.3.3 RSPM
cubature 2.1.4 RSPM
curl 6.2.3 CRAN (R 4.4.1)
CVXR 1.0-15 RSPM
dagitty 0.3-4 RSPM
data.table 1.17.4 RSPM
devtools 2.4.5 RSPM
DHARMa 0.4.7 RSPM
digest 0.6.37 RSPM
distributional 0.5.0 RSPM
doParallel 1.0.17 RSPM
doRNG 1.8.6.2 RSPM
doSNOW 1.0.20 RSPM
dplyr 1.1.4 RSPM
dreamerr 1.5.0 RSPM
ECOSolveR 0.5.5 RSPM
ellipsis 0.3.2 RSPM
emmeans 1.11.1 RSPM
estimability 1.5.1 RSPM
evaluate 1.0.3 RSPM
farver 2.1.2 RSPM
fastDummies 1.7.5 RSPM
fastmap 1.2.0 CRAN (R 4.4.3)
fect 1.0.0 RSPM
fixest 0.12.1 RSPM
forcats 1.0.0 RSPM
foreach 1.5.2 RSPM
forecast 8.24.0 RSPM
Formula 1.2-5 RSPM
fracdiff 1.5-3 RSPM
fs 1.6.6 RSPM
future 1.49.0 RSPM
future.apply 1.11.3 RSPM
geeM 0.10.1 RSPM
geepack 1.3.12 RSPM
generics 0.1.4 RSPM
geosphere 1.5-20 RSPM
GGally 2.2.1 RSPM
ggdag 0.2.13 RSPM
ggplot2 3.5.2 RSPM
ggstats 0.9.0 RSPM
glmmTMB 1.1.11 RSPM
glmnet 4.1-9 RSPM
glmx 0.2-1 RSPM
globals 0.18.0 RSPM
glue 1.8.0 RSPM
gmp 0.7-5 RSPM
gnm 1.1-5 RSPM
gridExtra 2.3 RSPM
gtable 0.3.6 RSPM
gtools 3.9.5 RSPM
hms 1.1.3 CRAN (R 4.4.3)
htmltools 0.5.8.1 RSPM
htmlwidgets 1.6.4 RSPM
httpuv 1.6.16 RSPM
httr2 1.1.2 RSPM
igraph 2.1.4 RSPM
iterators 1.0.14 RSPM
jsonlite 2.0.0 CRAN (R 4.4.3)
kableExtra 1.4.0 RSPM
kernlab 0.9-33 RSPM
knitr 1.50 RSPM
labeling 0.4.3 RSPM
later 1.4.2 RSPM
lattice 0.22-6 CRAN (R 4.4.1)
lfe 3.1.1 RSPM
lifecycle 1.0.4 RSPM
listenv 0.9.1 RSPM
lme4 1.1-37 RSPM
lmtest 0.9-40 RSPM
loo 2.8.0 RSPM
lubridate 1.9.4 RSPM
magrittr 2.0.3 RSPM
MASS 7.3-60.2 CRAN (R 4.4.1)
MatchIt 4.7.2 RSPM
Matrix 1.7-0 CRAN (R 4.4.1)
MatrixModels 0.5-4 RSPM
matrixStats 1.5.0 RSPM
maxLik 1.5-2.1 RSPM
memoise 2.0.1 CRAN (R 4.4.3)
mgcv 1.9-1 CRAN (R 4.4.1)
mime 0.13 CRAN (R 4.4.3)
miniUI 0.1.2 RSPM
minqa 1.2.8 RSPM
miscTools 0.6-28 RSPM
mvtnorm 1.3-3 RSPM
nixtlar 0.6.2 RSPM
nlme 3.1-164 CRAN (R 4.4.1)
nloptr 2.2.1 RSPM
nnet 7.3-19 CRAN (R 4.4.1)
np 0.60-18 RSPM
numDeriv 2016.8-1.1 RSPM
openxlsx 4.2.8 RSPM
optimx 2025-4.9 RSPM
PanelMatch 3.1.1 RSPM
parallelly 1.44.0 RSPM
pillar 1.10.2 RSPM
pkgbuild 1.4.8 RSPM
pkgconfig 2.0.3 RSPM
pkgload 1.4.0 RSPM
plm 2.6-6 RSPM
plyr 1.8.9 RSPM
posterior 1.6.1 RSPM
pracma 2.4.4 CRAN (R 4.4.1)
processx 3.8.6 RSPM
profvis 0.4.0 RSPM
promises 1.3.2 RSPM
ps 1.9.1 RSPM
purrr 1.0.4 RSPM
Qtools 1.5.9 RSPM
quadprog 1.5-8 RSPM
quantdr 1.2.2 RSPM
quantmod 0.4.28 RSPM
quantreg 6.1 RSPM
quarto 1.4.4 RSPM
qvcalc 1.0.4 RSPM
R6 2.6.1 RSPM
rappdirs 0.3.3 CRAN (R 4.4.3)
rbibutils 2.3 RSPM
RColorBrewer 1.1-3 RSPM
Rcpp 1.0.14 RSPM
RcppParallel 5.1.10 RSPM
Rdpack 2.6.4 RSPM
readr 2.1.5 CRAN (R 4.4.3)
reformulas 0.4.1 RSPM
relimp 1.0-5 RSPM
remotes 2.5.0 RSPM
renv 1.1.2 CRAN (R 4.4.1)
reshape2 1.4.4 RSPM
rgenoud 5.9-0.11 RSPM
rio 1.2.3 RSPM
rlang 1.1.6 RSPM
rmarkdown 2.29 RSPM
Rmpfr 1.1-0 RSPM
rngtools 1.5.2 RSPM
rstantools 2.4.0 RSPM
rstudioapi 0.17.1 RSPM
sandwich 3.1-1 RSPM
scales 1.4.0 RSPM
scpi 3.0.0 RSPM
sessioninfo 1.2.3 RSPM
shape 1.4.6.1 RSPM
shiny 1.10.0 RSPM
snow 0.4-4 RSPM
sp 2.2-0 RSPM
SparseM 1.84-2 RSPM
stringi 1.8.7 RSPM
stringmagic 1.2.0 RSPM
stringr 1.5.1 RSPM
survival 3.6-4 CRAN (R 4.4.1)
svglite 2.2.1 RSPM
Synth 1.1-8 RSPM
systemfonts 1.2.3 RSPM
tensorA 0.36.2.1 RSPM
textshaping 1.0.1 RSPM
tibble 3.2.1 RSPM
tidygraph 1.3.1 RSPM
tidyr 1.3.1 RSPM
tidyselect 1.2.1 RSPM
tidyverse 2.0.0 RSPM
timechange 0.3.0 RSPM
timeDate 4041.110 RSPM
TMB 1.9.17 RSPM
tseries 0.10-58 RSPM
TTR 0.24.4 RSPM
tzdb 0.5.0 CRAN (R 4.4.3)
urca 1.3-4 RSPM
urlchecker 1.0.1 RSPM
usethis 3.1.0 RSPM
V8 6.0.3 RSPM
vctrs 0.6.5 RSPM
viridisLite 0.4.2 RSPM
weathercan 0.7.3.9000 https://ropensci.r-universe.dev (R 4.4.3)
withr 3.0.2 RSPM
xfun 0.52 RSPM
xml2 1.3.8 CRAN (R 4.4.3)
xtable 1.8-4 RSPM
xts 0.14.1 RSPM
yaml 2.3.10 RSPM
zip 2.3.3 RSPM
zoo 1.8-14 RSPM
Code
reticulate::py_list_packages()%>% 
 knitr::kable(caption = "Python packages", format = "html",
      col.names = c("Package", "Version", "Requirement"),
    row.names = FALSE,
      align = c("c", "l", "r", "r"))%>% 
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12)|>
  kableExtra::scroll_box(width = "100%", height = "375px")  
Python packages
Package Version Requirement
absl-py 2.1.0 absl-py==2.1.0
asttokens 2.4.1 asttokens==2.4.1
astunparse 1.6.3 astunparse==1.6.3
audioconverter 2.0.3 audioconverter==2.0.3
autograd 1.6.2 autograd==1.6.2
autograd-gamma 0.5.0 autograd-gamma==0.5.0
beautifulsoup4 4.12.3 beautifulsoup4==4.12.3
Brotli 1.1.0 Brotli==1.1.0
certifi 2023.11.17 certifi==2023.11.17
cffi 1.16.0 cffi==1.16.0
charset-normalizer 3.3.2 charset-normalizer==3.3.2
clarabel 0.9.0 clarabel==0.9.0
click 8.1.7 click==8.1.7
cloudpickle 3.0.0 cloudpickle==3.0.0
colorama 0.4.6 colorama==0.4.6
comm 0.2.1 comm==0.2.1
contourpy 1.2.0 contourpy==1.2.0
cvxopt 1.3.2 cvxopt==1.3.2
cvxpy 1.5.2 cvxpy==1.5.2
cycler 0.12.1 cycler==0.12.1
debugpy 1.8.0 debugpy==1.8.0
decorator 4.4.2 decorator==4.4.2
delete-chrome-history-py 0.1.8 delete-chrome-history-py==0.1.8
easyocr 1.7.1 easyocr==1.7.1
ecos 2.0.13 ecos==2.0.13
editdistance 0.8.1 editdistance==0.8.1
efficientnet 1.0.0 efficientnet==1.0.0
essential-generators 1.0 essential-generators==1.0
et-xmlfile 1.1.0 et-xmlfile==1.1.0
executing 2.0.1 executing==2.0.1
fancyimpute 0.7.0 fancyimpute==0.7.0
ffmpeg 1.4 ffmpeg==1.4
ffmpeg-python 0.2.0 ffmpeg-python==0.2.0
filedir 0.0.3 filedir==0.0.3
filelock 3.13.1 filelock==3.13.1
flatbuffers 24.3.25 flatbuffers==24.3.25
fonttools 4.47.2 fonttools==4.47.2
formulaic 1.0.1 formulaic==1.0.1
fsspec 2023.12.2 fsspec==2023.12.2
future 0.18.3 future==0.18.3
gast 0.6.0 gast==0.6.0
git-filter-repo 2.45.0 git-filter-repo==2.45.0
google-pasta 0.2.0 google-pasta==0.2.0
graphviz 0.20.3 graphviz==0.20.3
grpcio 1.65.4 grpcio==1.65.4
gTTS 2.5.1 gTTS==2.5.1
h5py 3.11.0 h5py==3.11.0
idna 3.6 idna==3.6
imageio 2.34.2 imageio==2.34.2
imageio-ffmpeg 0.5.1 imageio-ffmpeg==0.5.1
imgaug 0.4.0 imgaug==0.4.0
iniconfig 2.0.0 iniconfig==2.0.0
interface-meta 1.3.0 interface-meta==1.3.0
ipykernel 6.29.5 ipykernel==6.29.5
ipython 8.20.0 ipython==8.20.0
jedi 0.19.1 jedi==0.19.1
Jinja2 3.1.3 Jinja2==3.1.3
joblib 1.4.0 joblib==1.4.0
jupyter_client 8.6.0 jupyter_client==8.6.0
jupyter_core 5.7.1 jupyter_core==5.7.1
keras 3.4.1 keras==3.4.1
Keras-Applications 1.0.8 Keras-Applications==1.0.8
keras-ocr 0.9.3 keras-ocr==0.9.3
kiwisolver 1.4.5 kiwisolver==1.4.5
knnimpute 0.1.0 knnimpute==0.1.0
lazy_loader 0.4 lazy_loader==0.4
libclang 18.1.1 libclang==18.1.1
lifelines 0.28.0 lifelines==0.28.0
llvmlite 0.41.1 llvmlite==0.41.1
Markdown 3.6 Markdown==3.6
markdown-it-py 3.0.0 markdown-it-py==3.0.0
MarkupSafe 2.1.4 MarkupSafe==2.1.4
matplotlib 3.8.2 matplotlib==3.8.2
matplotlib-inline 0.1.6 matplotlib-inline==0.1.6
mdurl 0.1.2 mdurl==0.1.2
mido 1.3.3 mido==1.3.3
ml-dtypes 0.4.0 ml-dtypes==0.4.0
more-itertools 10.2.0 more-itertools==10.2.0
moviepy 1.0.3 moviepy==1.0.3
mpmath 1.3.0 mpmath==1.3.0
multipledispatch 1.0.0 multipledispatch==1.0.0
mutagen 1.47.0 mutagen==1.47.0
namex 0.0.8 namex==0.0.8
natsort 8.4.0 natsort==8.4.0
nest-asyncio 1.5.9 nest-asyncio==1.5.9
networkx 3.2.1 networkx==3.2.1
ninja 1.11.1.1 ninja==1.11.1.1
nose 1.3.7 nose==1.3.7
numba 0.58.1 numba==0.58.1
numexpr 2.10.0 numexpr==2.10.0
numpy 1.26.3 numpy==1.26.3
openai-whisper 20231117 openai-whisper==20231117
opencv-python 4.10.0.84 opencv-python==4.10.0.84
opencv-python-headless 4.10.0.84 opencv-python-headless==4.10.0.84
openpyxl 3.1.4 openpyxl==3.1.4
opt-einsum 3.3.0 opt-einsum==3.3.0
optree 0.12.1 optree==0.12.1
osqp 0.6.5 osqp==0.6.5
packaging 23.2 packaging==23.2
pandas 2.2.0 pandas==2.2.0
pandas-flavor 0.6.0 pandas-flavor==0.6.0
parso 0.8.3 parso==0.8.3
patsy 0.5.6 patsy==0.5.6
pillow 10.2.0 pillow==10.2.0
platformdirs 4.1.0 platformdirs==4.1.0
pluggy 1.5.0 pluggy==1.5.0
polars 1.9.0 polars==1.9.0
proglog 0.1.10 proglog==0.1.10
prompt-toolkit 3.0.43 prompt-toolkit==3.0.43
protobuf 4.25.4 protobuf==4.25.4
psutil 5.9.8 psutil==5.9.8
pure-eval 0.2.2 pure-eval==0.2.2
pyarrow 15.0.0 pyarrow==15.0.0
pyclipper 1.3.0.post5 pyclipper==1.3.0.post5
pycparser 2.22 pycparser==2.22
pycryptodomex 3.20.0 pycryptodomex==3.20.0
pydotplus 2.0.2 pydotplus==2.0.2
pydub 0.24.1 pydub==0.24.1
Pygments 2.17.2 Pygments==2.17.2
pyjanitor 0.26.0 pyjanitor==0.26.0
PyMuPDF 1.24.9 PyMuPDF==1.24.9
PyMuPDFb 1.24.9 PyMuPDFb==1.24.9
pyparsing 3.1.1 pyparsing==3.1.1
PyPDF2 3.0.1 PyPDF2==3.0.1
pyreadr 0.5.0 pyreadr==0.5.0
pytesseract 0.3.10 pytesseract==0.3.10
pytest 8.3.1 pytest==8.3.1
python-bidi 0.6.0 python-bidi==0.6.0
python-dateutil 2.8.2 python-dateutil==2.8.2
pytube 15.0.0 pytube==15.0.0
pytube3 9.6.4 pytube3==9.6.4
pytz 2023.3.post1 pytz==2023.3.post1
pywin32 306 pywin32==306
PyYAML 6.0.1 PyYAML==6.0.1
pyzmq 25.1.2 pyzmq==25.1.2
qdldl 0.1.7.post1 qdldl==0.1.7.post1
regex 2023.12.25 regex==2023.12.25
requests 2.32.3 requests==2.32.3
rich 13.7.1 rich==13.7.1
rpy2 3.5.16 rpy2==3.5.16
scikit-image 0.24.0 scikit-image==0.24.0
scikit-learn 1.3.2 scikit-learn==1.3.2
scikit-survival 0.22.2 scikit-survival==0.22.2
scipy 1.11.4 scipy==1.11.4
scs 3.2.6 scs==3.2.6
seaborn 0.13.2 seaborn==0.13.2
semantic-version 2.10.0 semantic-version==2.10.0
setuptools-rust 1.8.1 setuptools-rust==1.8.1
shapely 2.0.5 shapely==2.0.5
six 1.16.0 six==1.16.0
soupsieve 2.5 soupsieve==2.5
SpeechRecognition 3.10.1 SpeechRecognition==3.10.1
spyder-kernels 2.5.2 spyder-kernels==2.5.2
stack-data 0.6.3 stack-data==0.6.3
statsmodels 0.14.1 statsmodels==0.14.1
sympy 1.12 sympy==1.12
target 0.0.11 target==0.0.11
tensorboard 2.17.0 tensorboard==2.17.0
tensorboard-data-server 0.7.2 tensorboard-data-server==0.7.2
tensorflow 2.17.0 tensorflow==2.17.0
tensorflow-intel 2.17.0 tensorflow-intel==2.17.0
tensorflow-io-gcs-filesystem 0.31.0 tensorflow-io-gcs-filesystem==0.31.0
termcolor 2.4.0 termcolor==2.4.0
threadpoolctl 3.4.0 threadpoolctl==3.4.0
tifffile 2024.7.24 tifffile==2024.7.24
tiktoken 0.5.2 tiktoken==0.5.2
torch 2.4.0 torch==2.4.0
torchaudio 2.4.0 torchaudio==2.4.0
torchvision 0.19.0 torchvision==0.19.0
tornado 6.4 tornado==6.4
tqdm 4.66.1 tqdm==4.66.1
traitlets 5.14.1 traitlets==5.14.1
translator 0.0.9 translator==0.0.9
typing_extensions 4.9.0 typing_extensions==4.9.0
tzdata 2023.4 tzdata==2023.4
tzlocal 5.2 tzlocal==5.2
urllib3 2.1.0 urllib3==2.1.0
validators 0.33.0 validators==0.33.0
watchdog 3.0.0 watchdog==3.0.0
wcwidth 0.2.13 wcwidth==0.2.13
websockets 12.0 websockets==12.0
Werkzeug 3.0.3 Werkzeug==3.0.3
whisper 1.1.10 whisper==1.1.10
wrapt 1.16.0 wrapt==1.16.0
xarray 2024.1.1 xarray==2024.1.1
youtube-dl 2021.12.17 youtube-dl==2021.12.17
yt-dlp 2024.7.9 yt-dlp==2024.7.9

References

1.
Brodersen KH, Gallusser F, Koehler J, Remy N, Scott SL. Inferring causal impact using bayesian structural time-series models. The Annals of Applied Statistics. 2015;9(1). doi:10.1214/14-aoas788.